import numpy as np
import matplotlib.pyplot as plt
import pyproj
import dolfin
import folium
import femlium
from utils import gmsh_to_fenics
Auxiliary function to get a folium Map close to Lake Garda.
def get_garda_geo_map(boundary_icons=False):
# Add map close to Lake Garda
geo_map = folium.Map(location=[45.6389113, 10.7521368], zoom_start=10.3)
# Add markers
if boundary_icons:
location_markers = {
"Sarca": [45.87395405, 10.87087005],
"Mincio": [45.43259035, 10.7007715]
}
location_colors = {
"Sarca": "red",
"Mincio": "green"
}
for key in location_markers.keys():
folium.Marker(
location=location_markers[key],
tooltip=key,
icon=folium.Icon(color=location_colors[key])
).add_to(geo_map)
# Return folium map
return geo_map
get_garda_geo_map()
Read the mesh, the subdomain markers and the boundary markers from file with dolfin.
mesh, subdomains, boundaries = gmsh_to_fenics("data/garda.msh")
Plot the mesh using dolfin.plot.
fig = plt.figure(figsize=(12, 12))
dolfin.plot(mesh)
fig.gca().axis("equal")
(616658.03914899996, 647061.96085100004, 5029877.2258409997, 5085382.7741590003)
Define a pyproj Transformer to map between different reference systems, because the points read from file are stored a $(x, y)$ pairs in the EPSG32632 reference system, while the map produced by folium is based on (latitude, longitude) pairs in the EPSG4326 reference system.
transformer = pyproj.Transformer.from_crs("epsg:32632", "epsg:4326", always_xy=True)
We define a mesh plotter for meshes in dolfin format, which is implemented in femlium.DolfinPlotter.
dolfin_plotter = femlium.DolfinPlotter(transformer)
We use the dolfin_plotter to draw the mesh on top of the geographic map.
geo_map = get_garda_geo_map()
dolfin_plotter.add_mesh_to(geo_map, mesh)
geo_map
We may change the color and the weight of the line.
geo_map = get_garda_geo_map()
dolfin_plotter.add_mesh_to(geo_map, mesh, face_colors="red", face_weights=2)
geo_map
Furthermore, we may set the colors and the weights of the face representation to depend on the markers associated to each segment.
geo_map = get_garda_geo_map(boundary_icons=True)
face_colors = {
0: "gray",
1: "blue",
2: "red",
3: "green"
}
face_weights = {
0: 1,
1: 2,
2: 5,
3: 5
}
dolfin_plotter.add_mesh_to(
geo_map, mesh, face_mesh_function=boundaries, face_colors=face_colors, face_weights=face_weights)
geo_map
Cells can be colored as well, with a uniform color or depending on the cell markers. We start from a uniform color.
geo_map = get_garda_geo_map()
dolfin_plotter.add_mesh_to(geo_map, mesh, cell_colors="orange")
geo_map
We also show the case of colors being set from cell markers. There are two cell markers in this mesh, equal to 1 for the region close to the shoreline (colored in purple) and 2 for the rest of the domain (colored in yellow).
geo_map = get_garda_geo_map()
cell_colors = {
1: "purple",
2: "yellow"
}
dolfin_plotter.add_mesh_to(geo_map, mesh, cell_mesh_function=subdomains, cell_colors=cell_colors)
geo_map
Once can use colors associated to both cell and face markers on the same plot.
geo_map = get_garda_geo_map(boundary_icons=True)
dolfin_plotter.add_mesh_to(
geo_map, mesh,
cell_mesh_function=subdomains, face_mesh_function=boundaries,
cell_colors=cell_colors, face_colors=face_colors, face_weights=face_weights)
geo_map
In order to define a simple scalar field, we compute the centroid of the domain.
centroid = np.mean(mesh.coordinates(), axis=0)
We may plot the centroid on top of the mesh.
geo_map = get_garda_geo_map()
dolfin_plotter.add_mesh_to(geo_map, mesh)
folium.Marker(location=transformer.transform(*centroid)[::-1], tooltip="Centroid").add_to(geo_map)
geo_map
The scalar field is defined as $s(\rho, \theta) = \frac{\rho}{\sqrt{1 - 0.5 \cos^2 \theta}}$, and is interpolated on a $\mathbb{P}^2$ finite element space. Here $(\rho, \theta)$ are the polar coordinates centered at the centroid.
scalar_function_space = dolfin.FunctionSpace(mesh, "CG", 2)
Calling FFC just-in-time (JIT) compiler, this may take some time.
class ScalarField(dolfin.UserExpression):
def eval_cell(self, value, x, cell):
rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
value[0] = rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)
def value_shape(self):
return ()
scalar_field = dolfin.interpolate(ScalarField(), scalar_function_space)
We next show a filled contour plot with 15 levels using dolfin.plot.
fig = plt.figure(figsize=(12, 12))
trif = dolfin.plot(scalar_field, mode="contourf", levels=15, cmap="jet")
fig.colorbar(trif)
fig.gca().axis("equal")
(618040.03558999998, 645679.96441000002, 5032400.2053100001, 5082859.7946899999)
In order to plot a field on a geographic map, we use again the dolfin_plotter. We may plot the same filled contour plot on the geographic map.
geo_map = get_garda_geo_map()
dolfin_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contourf", levels=15, cmap="jet")
geo_map
Similarly, we can also use (unfilled) contour plots.
fig = plt.figure(figsize=(12, 12))
tri = dolfin.plot(scalar_field, mode="contour", levels=15, cmap="jet")
fig.colorbar(tri)
fig.gca().axis("equal")
(618040.03558999998, 645679.96441000002, 5032400.2053100001, 5082859.7946899999)
geo_map = get_garda_geo_map()
dolfin_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map
One may also combine mesh plots and solution plots.
geo_map = get_garda_geo_map()
dolfin_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
dolfin_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map
We next define a vector field $\mathbf{v}(\rho, \theta) = \begin{bmatrix}-\rho \sin \theta\\\rho \cos\theta \end{bmatrix}$.
vector_function_space = dolfin.VectorFunctionSpace(mesh, "CG", 2)
Calling FFC just-in-time (JIT) compiler, this may take some time.
class VectorField(dolfin.UserExpression):
def eval_cell(self, value, x, cell):
rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
value[0] = - rho * np.sin(theta)
value[1] = rho * np.cos(theta)
def value_shape(self):
return (2, )
vector_field = dolfin.interpolate(VectorField(), vector_function_space)
We may obtain contourf or contour plots of the magnitude of the vector field.
geo_map = get_garda_geo_map()
dolfin_plotter.add_vector_field_to(geo_map, vector_field, mode="contourf", levels=15, cmap="jet")
geo_map
geo_map = get_garda_geo_map()
dolfin_plotter.add_vector_field_to(geo_map, vector_field, mode="contour", levels=15, cmap="jet")
geo_map
Vector field can also be plotted using a quiver. We first see the quiver plot obtained with dolfin.plot.
fig = plt.figure(figsize=(12, 12))
quiv = dolfin.plot(vector_field, mode="glyphs", cmap="jet")
fig.colorbar(quiv)
fig.gca().axis("equal")
(616658.03914899996, 647061.96085100004, 5029877.2258409997, 5085382.7741590003)
A similar plot can rendered on top of the geographic map.
geo_map = get_garda_geo_map()
dolfin_plotter.add_vector_field_to(geo_map, vector_field, mode="quiver", scale=1e-1, cmap="jet")
geo_map